function outparam = fitplot(data, fit_func, clr)

global f1 f2 fbeat RabiTime gSense gFitStarter

if length(data) == 3
    xdata = data{1};
    ydata = data{2};
    weights = data{3};
    xmin = min(xdata);
    xmax = max(xdata);
elseif length(data) == 2
    xdata = data{1};
    ydata = data{2};
    weights = ones(length(xdata),1);
    xmin = min(xdata);
    xmax = max(xdata);
elseif length(data) == 1
    xdata = 1:size(data{1},2);
    ydata = data{1};
    weights = ones(length(xdata),1);
    xmin = min(xdata);
    xmax = max(xdata);
elseif length(data) == 4
    xdata = data{1};
    ydata = data{2};
    weights = data{3};
    xmin = min(xdata);
    xmax = max(xdata);
end

if size(xdata,1) == 1
    xdata = xdata';
end
if size(ydata,1) == 1
    ydata = ydata';
end

switch fit_func
    case 'strexp0.5'
        FitType = fittype('a*exp(-sqrt(x/b))');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [ydata(1) xdata(end)/2], 'lower', [-inf 0],'weights', 1./weights, 'robust', 'BiSquare');
    case 'logstrexp'
        FitType = fittype('log(a*exp(-(x/b).^c))');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [ydata(1) 0.38 0.5], 'lower', [0 0 0], 'upper', [inf inf inf],'weights', 1./weights, 'robust', 'BiSquare');
    case 'strexp'
        FitType = fittype('a*exp(-(x/b).^c)');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [ydata(1) xdata(end)/2 0.5], 'lower', [ydata(1)*0.95 0 0], 'upper', [ydata(1)*1.05 inf inf],'weights', 1./weights, 'robust', 'BiSquare');
    case 'strexp_rate'
        FitType = fittype('a*exp(-(b*x).^c)');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [ydata(1) xdata(end)/2 0.5], 'lower', [0 0 0], 'upper', [inf inf inf],'weights', 1./weights, 'robust', 'BiSquare');
    case 'strexp_offset'
        FitType = fittype('a*exp(-(x/b).^c) + d');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [ydata(1) xdata(end)/2 1 0], 'lower', [0 0 0 0], 'upper', [inf inf inf inf],'weights', 1./weights, 'robust', 'BiSquare');
    case 'strexp_modulated'
        if isempty(gFitStarter)
           gFitStarter = rand; 
        end
        FitType = fittype('a*exp(-(x/b).^c).*(e*cos(2*pi*x/d)+(1-e))');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [ydata(1) xdata(end)/2 1 gFitStarter 1], 'lower', [ydata(1)*0.99 0 0 0 0], 'upper', [ydata(1)*1.01 inf 1 inf 1],'weights', 1./weights, 'robust', 'BiSquare');
    case 'strexp_modulated_offset'
        if isempty(gFitStarter)
            gFitStarter = rand;
        end
        FitType = fittype('a*exp(-(x/b).^c).*(e*cos(2*pi*x/d)+(1-e)) + f');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [ydata(1)-gFitStarter xdata(end)/2 1 10 1 gFitStarter], 'lower', [(ydata(1)-gFitStarter)*0.99 0 0 0 0 gFitStarter*0.999], 'upper', [(ydata(1)-gFitStarter)*1.01 inf 1 inf 1 gFitStarter*1.001],'weights', 1./weights, 'robust', 'BiSquare');
    case 'superexp'
        FitType = fittype('a*exp(-(x/b).^c)');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [ydata(1) xdata(end)/2 2], 'lower', [0 0 1], 'upper', [inf inf inf],'weights', 1./weights);
    case 'superexp_rate'
        FitType = fittype('a*exp(-(b*x).^c)');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [ydata(1) xdata(end)/2 2], 'lower', [0 0 1], 'upper', [inf inf inf],'weights', 1./weights, 'robust', 'BiSquare');
    case 'simpleexp'
        FitType = fittype('a*exp(-(x/b))');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [ydata(1) xdata(end)/2], 'lower', [-inf 0], 'upper', [inf inf],'weights', 1./weights);
    case 'simpleexp_offset'
        FitType = fittype('a*exp(-(x/b)) + c');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [ydata(1) xdata(end)/2 0], 'lower', [-inf -inf -inf], 'upper', [inf inf inf],'weights', 1./weights);
    case 'simpleexp_rate'
        FitType = fittype('a*exp(-(b*x))');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [ydata(1) 1/(xdata(end)/2)], 'lower', [0 0], 'upper', [inf inf], 'weights', 1./weights, 'robust', 'BiSquare');
    case 'Arrhenius'
        FitType = fittype('a*exp(b*x)');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [ydata(1) 1/(xdata(end)/2)], 'lower', [0 0], 'upper', [inf inf], 'weights', 1./weights, 'robust', 'BiSquare');
    case 'doubleexp'
        FitType = fittype('log(a*exp(-(x/b).^2) + c*exp(-(x/d)))');
        FitRes = fit(xdata, log(ydata), FitType, 'startpoint', [ydata(1)*0.9 xdata(end)*0.1 ydata(1)*0.1 xdata(end)*0.9], 'lower', [0 0.01 0 0.01], 'upper', [inf inf inf inf],'weights', 1./weights);
    case 'doubleexp1'
        FitType = fittype('a*exp(-(x/b)) + c*exp(-(x/d))');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [ydata(1)*0.9 xdata(end)*0.1 ydata(1)*0.1 xdata(end)*0.9], 'lower', [0 0.01 -inf 0.01], 'upper', [inf inf inf inf],'weights', 1./weights);
    case 'doubleexp_rate'
        FitType = fittype('a*exp(-((b+d)*x)) + c*exp(-(d*x))');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [ydata(1)*0.9 1/(xdata(end)*1) ydata(1)*0.1 1/(xdata(end)*1)], 'lower', [0 0 0 0], 'upper', [inf inf inf inf],'weights', 1./weights);
    case 'linear'
        FitType = fittype('a*x+b');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [1 -10], 'weights', 1./weights);
    case 'linear_mirror'
        FitType = fittype('a*abs(x-b)');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [1e3 0], 'weights', 1./weights);
    case 'linear0'
        FitType = fittype('a*x');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [1], 'weights', 1./weights, 'robust', 'BiSquare');
    case 'linearkink'
        FitType = fittype('a*(x-c)*heaviside(x-c)+b');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [10 0.15 0], 'lower', [-inf -inf 0], 'upper', [inf inf 20], 'weights', 1./weights, 'robust', 'BiSquare');
    case 'two_linears'
        FitType = fittype('a*x+b');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [1 -10]);
    case '1/x'
        FitType = fittype('a./x');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [1]);
    case '1/x2'
        FitType = fittype('a./x.^2');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [10], 'robust', 'BiSquare');
    case '1/xsqrt'
        FitType = fittype('a./x.^0.5');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [10], 'robust', 'BiSquare');
    case 'power'
        FitType = fittype('a*x.^b');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [10 0.5], 'lower', [0 0], 'upper', [inf 1], 'robust', 'BiSquare');
    case 'U-shaped'
        FitType = fittype('a*abs((x-b)).^c + d');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [3 0 2 0.12], 'lower', [0 -inf 0 0], 'upper', [inf inf inf inf], 'robust', 'BiSquare');
    case 'quadratic'
        FitType = fittype('a*(x-b).^2 + c');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [1e5 mean(xdata) min(ydata)], 'robust', 'BiSquare');
    case 'quadratic-asym'
        FitType = fittype('(x-b).^2.*(a1*heaviside(x-b)+a2*heaviside(b-x)) + c');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [1/2 1/2 0 0], 'lower', [0 0 0.999 -inf], 'upper', [inf inf 1.001 inf], 'robust', 'BiSquare');
    case 'quadratic-1/2'
        FitType = fittype('1/2*pi^2*x.^2 + a');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [0], 'lower', [0], 'upper', [inf], 'robust', 'BiSquare');
    case 'polyn'
        FitType = fittype('a*pi^2*x.^b + c');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [1 3 0], 'lower', [0 2 0], 'robust', 'BiSquare');
    case 'exp_positive'
        FitType = fittype('a*exp(b*x)');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [10 5]);
    case 'cubic'
        FitType = fittype('a*x + b*x.^2 + c*x.^3');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [10 1 1]);
    case 'Gaussian'
        FitType = fittype('a*exp(-0.5*((x-b)/c).^2) + d');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [max(ydata) (xdata(1)+xdata(end))/3 abs(xdata(end)-xdata(1))*0.2 0], 'lower', [-inf -inf -inf -inf], 'upper', [inf inf inf inf],'weights', 1./weights);
    case 'InvGaussian'
        FitType = fittype('d - a*exp(-0.5*((x-b)/c).^2)');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [max(ydata) (xdata(1)+xdata(end))/2 abs(xdata(end)-xdata(1))*0.2 1], 'lower', [-inf -inf 0 -inf], 'upper', [inf inf inf inf],'weights', 1./weights);
    case 'NormInvGaussian'
        FitType = fittype('1 - a*exp(-0.5*((x-b)/c).^2)');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [max(ydata) (xdata(1)+xdata(end))/2 abs(xdata(end)-xdata(1))*0.2], 'lower', [-inf -inf -inf], 'upper', [inf inf inf],'weights', 1./weights);
    case 'SuperGaussian'
        FitType = fittype('a*exp(-0.5*(x/b).^c) + d');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [max(ydata) abs(xdata(end)-xdata(1))*0.2 2 0], 'lower', [0 0 0 0], 'upper', [inf inf inf 0.05],'weights', 1./weights);
    case 'AsymSuperGaussian'
        FitType = fittype('a*(heaviside(x-d)*exp(-0.5*(abs((x-d)/b1)).^c) + heaviside(d-x)*exp(-0.5*(abs((x-d)/b2)).^c))');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [0.4 0.01 0.01 1 1], 'lower', [0 0 0 1 -0.05], 'upper', [1 1 1 inf 0.05], 'weights', 1./weights);
    case 'PhaseTransition'
        FitType = fittype('a*(heaviside(x-b1).*heaviside(b2-x)) .* ((x-b1).^(c1) .* (b2-x).^(c2))');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [0.5 0.98 1.02 0.3 0.3], 'lower', [1 0.7 1 0 0], 'upper', [1 1 1.3 10 10], 'weights', 1./weights);
    case '3T_ansatz'
        FitType = fittype(['c*heaviside(b-x) + (a*(x-b).^2 + c).*heaviside(x-b)']);
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [1 0.05 0], 'lower', [0 0 0], 'upper', [inf 0.2 inf], 'weights', 1./weights, 'robust', 'BiSquare');
    case 'critical_DTC'
        FitType = fittype(['a*(x*pi).^2 .*exp(-b*(x*pi).^(-2)) + c']);
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [0.5 0.1 0.01], 'lower', [0 0 0], 'upper', [inf inf inf], 'weights', 1./weights, 'robust', 'BiSquare');
    case '3T_Phase_int'
        FitType = fittype(['a/(6*c*pi^1/2)*(exp(-b./(x+c).^2/pi^2).*(sqrt(pi).*(x+c).*(-2*b+pi^2*(c+x).^2) - 2*b^(3/2).*exp(b./(x+c).^2/pi^2).*erf(sqrt(b)/pi./(x+c)))+exp(-b./(c-x).^2/pi^2).*(sqrt(pi).*(c-x).*(-2*b+pi^2*(c-x).^2) - 2*b^(3/2).*exp(b./(x-c).^2/pi^2).*erf(sqrt(b)/pi./(c-x))))']);
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [1 0.1 0.005], 'lower', [0 0 0.001], 'upper', [inf inf 0.008], 'weights', 1./weights);
    case '3T_Phase_tau'
        FitType = fittype(['a*(x*0.06*pi).^(1).*exp(-b*(0.06*pi)^(-2)*(x).^(-1)) + c']);
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [1 80 0.01], 'lower', [0 10 0], 'upper', [1 300 0.1], 'weights', 1./weights);
    case 'Lorentzian'
        FitType = fittype('a./(((x-b)/c).^2 + 1) + d');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [max(ydata) (xdata(1)+xdata(end))/2 abs(xdata(end)-xdata(1))*0.2 0], 'lower', [-inf -inf -inf -inf], 'upper', [inf inf inf inf],'weights', 1./weights);
    case '2Lorentzian'
        FitType = fittype('a./(((x-b)/c).^2 + 1) + e./(((x-f)/g).^2 + 1) + d');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [-max(ydata)*0.0015 (xdata(1)+xdata(end))/2.16 abs(xdata(end)-xdata(1))/5.8 1 -max(ydata)*0.002 (xdata(1)+xdata(end))/1.9 abs(xdata(end)-xdata(1))/5.8], 'lower', [-inf min(xdata) -inf -inf -inf min(xdata) -inf], 'upper', [inf max(xdata) inf inf inf max(xdata) inf],'weights', 1./weights);
    case 'InvLorentzian'
        FitType = fittype('d - a./(((x-b)/c).^2 + 1)');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [max(ydata) (xdata(1)+xdata(end))/2 abs(xdata(end)-xdata(1))*0.2 max(ydata)-min(ydata)], 'lower', [0 -inf 0 -inf], 'upper', [inf inf inf inf],'weights', 1./weights);   
    case 'NormInvLorentzian'
        FitType = fittype('1 - a./(((x-b)/c).^2 + 1)');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [max(ydata) (xdata(1)+xdata(end))/2 abs(xdata(end)-xdata(1))*0.2], 'lower', [0 -inf 0], 'upper', [inf inf inf],'weights', 1./weights); 
    case 'NormInvLorentzianSplitting'
        FitType = fittype('1 - heaviside(x-b).*a./(((x-(b+d))/c).^2 + 1) - heaviside(b-x).*a./(((x-(b-d))/c).^2 + 1)');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [max(ydata) (xdata(1)+xdata(end))/2 abs(xdata(end)-xdata(1))*0.2 (xdata(1)+xdata(end))/1000], 'lower', [0 -inf 0 0], 'upper', [inf inf inf (xdata(1)+xdata(end))/2*0.01],'weights', 1./weights);
    case 'NormInvTwoLorentzians'
        FitType = fittype('1 - a1./(((x-b1)/c1).^2 + 1) - a2./(((x-b2)/c2).^2 + 1)');
        starter = [0.3,0.1,0.746,0.73,0.005,0.005];
        FitRes = fit(xdata, ydata, FitType, 'startpoint', starter,...
            'lower', [0 0 -inf -inf 0 0], 'upper', [inf inf inf inf inf 0.01],'weights', 1./weights);   
    case 'InvTwoLorentzians'
        FitType = fittype('d - a1./(((x-b1)/c1).^2 + 1) - a2./(((x-b2)/c2).^2 + 1)');
        starter = [0.3,0.1,0.749,0.737,0.005,0.005,1];
        lb = [0 0 0 0 0 0 0.5];%[0 0 0.748 0.736 0 0 0.5];
        ub = [inf inf inf inf inf 0.01 1.1];%[inf inf 0.75 0.738 inf 0.01 1.1];
        FitRes = fit(xdata, ydata, FitType, 'startpoint', starter,...
            'lower', lb, 'upper', ub,'weights', 1./weights);   
    case 'TwoLorentzians'
        dw = abs(xdata(1)-xdata(end));
        FitType = fittype('a./(((x-b)/c).^2 + 1) + d./(((x-f)/e).^2 + 1)');
        FitRes = fit(xdata, ydata, FitType,...
            'startpoint', [max(ydata) (xdata(1)+xdata(end))/2 0.2*abs(xdata(end)-xdata(1)) ... 
                           max(ydata) 0.2*abs(xdata(end)-xdata(1)) (xdata(1)+xdata(end))/2],...
            'lower', [0             -inf    0.01*dw      0            0.1*dw    -inf],...
            'upper', [max(ydata)    inf     inf          max(ydata)   inf       inf],...
            'weights', 1./weights);
    case 'ThreeLorentzians'
        dw = abs(xdata(1)-xdata(end));
        FitType = fittype('a./(((x-b)/c).^2 + 1) + d./(((x-0.5-e)/f).^2 + 1) + d./(((x-0.5+e)/f).^2 + 1)');
        FitRes = fit(xdata, ydata, FitType,...
             'startpoint', [max(ydata)      (xdata(1)+xdata(end))/2 abs(xdata(end)-xdata(1))*0.2 ... 
                            0.5*max(ydata)  0.025                   0.02],...
            'lower', [0                 0   0          0          0     0],...
            'upper', [max(ydata)        inf 0.01*dw    max(ydata) inf   dw],...
            'weights', 1./weights);
    case 'NormInvSinc'
        FitType = fittype('1 - a*sinc(x-b)');
    case 'OscDecay'     % oscillation and decay
        FitType = fittype('a*exp(-(x/b).^1)*(1+e*cos(2*pi*d*x))/(1+e)');
        FitRes = fit(xdata, ydata, FitType, 'startpoint',starter,...
            'lower', [0 0 0 0], 'upper', [inf inf inf inf],...
            'weights', 1./weights, 'robust', 'BiSquare'); 
    case 'OscDecayPhase'     % oscillation and decay
        FitType = fittype('a*exp(-(x/b).^1)*(1+e*cos(2*pi*d*x+f))/(1+e)');
        FitRes = fit(xdata, ydata, FitType, 'startpoint',...
            [max(abs(ydata)),(xdata(end)-xdata(1))/1,0.7/(xdata(end)-xdata(1)),2,0],...
            'lower', [-inf 0 0 0 -pi], 'upper', [inf inf inf inf pi],...
            'weights', 1./weights, 'robust', 'BiSquare'); 
    case 'OscDecayStr'     % oscillation and decay
        FitType = fittype('a*exp(-(x/b).^c)*(1+e*cos(2*pi*d*x))/(1+e)');
        FitRes = fit(xdata, ydata, FitType, 'startpoint',...
            [max(ydata),(xdata(end)-xdata(1))/2,1,2/(xdata(end)-xdata(1)),1000],...
            'lower', [0 0 0 0 0], 'upper', [0.1 inf 5 inf inf],...
            'weights', 1./weights, 'robust', 'BiSquare'); 
    case 'Rabi'
        FitType = fittype('1 - a/2*(1-exp(-x/c).*cos(2*pi*b*x))');
        [~,loc] = findpeaks(ydata);
        if isempty(loc)
            loc = xdata(end);
        else
            loc = loc(1);
        end
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [(1-min(ydata))*2 1/RabiTime 500], 'lower', [0 0 0], 'upper', [3 inf inf], 'weights', 1./weights, 'robust', 'BiSquare'); 
    case 'QST'
        FitType = fittype('(d+a)/2 + (d-a)/2*cos(2*pi*c*x) + b*sin(2*pi*c*x)');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [0.5 0.1 1/14 0.5], 'lower', [0 -inf 1/16 0], 'upper', [1 inf 1/13 1], 'weights', 1./weights, 'robust', 'BiSquare'); 
    case 'RabiShift'
        FitType = fittype('1 - d + a*exp(-x/c).*(1-cos(2*pi*b*x+e))');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [0.05 1/RabiTime 500 0 0.05], 'lower', [0 0 10 0 -pi], 'upper', [1 inf inf 1 pi], 'weights', 1./weights, 'robust', 'BiSquare');
    case 'DoubleRabi'
        FitType = fittype('1 - a1/2*(1-exp(-x/c).*cos(2*pi*b1*x)) - a2/2*(1-exp(-x/c).*cos(2*pi*b2*x))');
        [~,loc] = findpeaks(ydata);
        loc = loc(1);
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [0.05 0.05 80e-3 80e-3 500], 'lower', [0 0 60e-3 60e-3 0], 'upper', [0.1 0.1 100e-3 100e-3 inf], 'weights', 1./weights);
    case 'opt_lifetime'
        FitType = fittype('a * exp(-((b-x)/c).^2)*heaviside(b-x) + a * (f*exp(-((x-b)/d)) + (1-f)*exp(-((x-b)/e))) * heaviside(x-b)');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [14 40 1 10 50 0.5], 'lower', [0 0 0 0 0 0], 'upper', [inf inf inf inf inf 1], 'weights', 1./weights);
    case 'powerlaw'
        FitType = fittype('a*x.^(-b)');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [1 0.3], 'lower', [0 0], 'upper', [inf inf], 'weights', 1./weights);
    case 'powerlawext'
        FitType = fittype('a * heaviside(b-x) + (c * (x-b).^d + a) .* heaviside(x-b)');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [0.1 2 0.01 0.75], 'lower', [0 0 0 0.5], 'upper', [inf inf inf 1], 'weights', 1./weights);
    case 'DTCdecay'
        FitType = fittype('log(a1*exp(-b1*x) + a2*exp(-b2*x) + a3*exp(-b1*x)*cos(2*pi*c*x))');
        FitRes = fit(xdata, log(ydata), FitType, 'startpoint', [50 5 20 1/15 1/40 fbeat], 'lower', [0 0 0 0 0 fbeat*0.9], 'upper', [inf inf inf 1 1 fbeat*1.1], 'weights', 1./weights, 'robust', 'BiSquare');        
    case 'ThreePeakDTC'
        FitType = fittype('abs(a1/((x-b1)/c1 - 1i)).^2 + abs(a2/((x-1/3)/c2 - 1i)).^2 + abs(a3/((x-b2)/c3 - 1i)).^2'); 
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [1 2 3 f1 f2 0.1 0.1 0.1], 'lower', [0 0 0 f1*0.999 f2*0.999 0 0 0], 'upper', [10 10 10 f1*1.001 f2*1.001 0.1 0.1 0.1], 'weights', 1./weights);
    case 'asymThreePeakDTC'
        fitfunc1 = 'abs(heaviside(x-b1) * a1/((x-b1)/c1 - 1i) + heaviside(b1-x) * a1/((x-b1)/d1 - 1i) + a2/((x-1/3)/c2 - 1i) + heaviside(x-b2) * a3/((x-b2)/c3 - 1i) + heaviside(b2-x) * a3/((x-b2)/d2 - 1i)).^2';
        FitType = fittype(fitfunc1);
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [1 10 3 f1 f2 0.1 0.1 0.1 0.1 0.1], 'lower', [0 0 0 f1*0.999 f2*0.999 0.005 0.005 0.005 0.005 0.005], 'upper', [10 10 10 f1*1.001 f2*1.001 0.1 0.1 0.1 0.1 0.1], 'weights', 1./weights);
    case 'sin1'
        FitType = fittype(fit_func);
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [max(abs(ydata)) 1/(max(xdata)-min(xdata)) 0]);
    case 'sinfreq'
        FitType = fittype(['a1*sin(',num2str(f1),'*x+c1)']);
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [max(abs(ydata)) 1]);
    case 'sinfreqshift'
        FitType = fittype(['a1*sin(',num2str(f1),'*x+c1)+d1']);
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [max(abs(ydata)) 1 0],'lower',[0,-inf,-inf]);
    case 'cos'
        FitType = fittype(['a1*cos(b1*x)']);
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [max(abs(ydata)) 5/(max(xdata)-min(xdata))]);
    case 'cosshift'
        FitType = fittype(['a1*cos(b1*x)+d1']);
        FitRes = fit(xdata, ydata, FitType, 'startpoint',...
            [(max(ydata)-min(ydata))/2,2/(max(xdata)-min(xdata)),(max(abs(ydata))+min(abs(ydata)))/2]);
    case 'cosshiftphase'
        FitType = fittype(['a1*cos(b1*x+c1)+d1']);
        FitRes = fit(xdata, ydata, FitType, 'startpoint',...
            [(max(ydata)-min(ydata))/2,2*pi/(max(xdata)-min(xdata)),0,(max(ydata)+min(ydata))/2],...
            'lower',[0.0,0,-pi,-inf],'upper',[inf,inf,pi,inf]);
    case 'normcosshift'
        FitType = fittype(['a1*cos(b1*x)+1-a1']);
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [(1-min(abs(ydata)))/2 5/(max(xdata)-min(xdata))]);
    case 'coscosfreq'
        FitType = fittype(['cos(b1*cos(',num2str(f1),'*x-c1))']);
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [5/(max(xdata)-min(xdata)) 0],...
            'lower',[-inf,-pi],'upper',[inf,pi]);
    case 'coscosfreqamp'
        FitType = fittype(['a1*cos(b1*cos(',num2str(f1),'*x-c1))+1-a1']);
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [(1-min(ydata))/2 1*pi/(max(ydata)-min(ydata)) 0],...
            'lower',[0,-inf,-pi],'upper',[1,inf,pi]);
    case 'twosinfreq'
        FitType = fittype(['a1*sin(',num2str(f1(1)),'*x+c1)+a2*sin(',num2str(f1(2)),'*x+c2)+d']);
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [max(abs(ydata))/2 0 max(abs(ydata))/2 0 0]);
    case 'ThirdOrder'
        FitType = fittype('a*x + b*x^2 + c*x^3');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [0.01 0 0], 'lower', [-inf -inf -inf], 'upper', [inf inf 300]);
    case 'ThirdOrderOdd'
        FitType = fittype('a*x + c*x^3');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [0.01 0], 'lower', [-inf -inf], 'upper', [inf 300]);
    case 'SatFit'
        FitType = fittype('a*x/(b+x)');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [0.01,1]);
    case 'cos_with_offset'
        FitType = fittype('a*cos(4*pi*x+b)+c'); 
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [(max(ydata)-min(ydata))/2 0 0], 'lower', [0 0 -inf], 'upper', [max(ydata)-min(ydata) 2*pi inf], 'weights', 1./weights, 'robust', 'BiSquare');
    case 'coscos_with_offset'
        FitType = fittype('a*cos(b*cos(2*pi*x+c))+d');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [(max(ydata)-min(ydata)) 1 pi 0], 'lower', [0 0 -2*pi -inf], 'upper', [10*(max(ydata)-min(ydata)) inf 2*pi inf], 'weights', 1./weights, 'robust', 'BiSquare');
    case 'sinewavedecay'
        FitType = fittype('a*exp(-((x+d)/c).^2).*sin(2*pi/b*(x+d))+e');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [gSense(1)*(max(ydata)-min(ydata)) gSense(2) gSense(3) 0 gSense(4)], 'lower', [-inf 0 0 -inf -inf], 'upper', [inf inf inf inf inf], 'weights', 1./weights, 'robust', 'BiSquare');
    case 'sinewavedecayperfect'
        FitType = fittype('a*exp(-(x/c).^2).*sin(2*pi/b*x)+e');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [gSense(1)*(max(ydata)-min(ydata)) gSense(2) gSense(3) gSense(4)], 'lower', [-inf 0 0 -inf], 'upper', [inf inf inf inf], 'weights', 1./weights, 'robust', 'BiSquare');
    case 'sinewave'
        FitType = fittype('a*sin(2*pi/b*x)');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [gSense(1)*(max(ydata)-min(ydata))/2 gSense(2)], 'lower', [-inf 0], 'upper', [inf inf], 'weights', 1./weights, 'robust', 'BiSquare');
    case 'sinewaveoffset'
        FitType = fittype('a*sin(2*pi/b*x+c)+d');
        FitRes = fit(xdata, ydata, FitType, 'startpoint', [gSense(1)*(max(ydata)-min(ydata))/2 1e3*gSense(2) 0 0], 'lower', [-inf 0 -pi -inf], 'upper', [inf inf pi inf], 'weights', 1./weights, 'robust', 'BiSquare');
end

% set plot properties manually
if strcmp(fit_func, 'powerlaw')
    xfit = linspace(xmin, xmax, 1000);
    linestyle = '-';
    linewidth = 0.5;
elseif strcmp(fit_func, 'cubic')% || strcmp(fit_func, 'quadratic')
    xfit = linspace(0, xmax, 1000);
    linestyle = '-';
    linewidth = 0.5;
elseif  strcmp(fit_func, 'quadratic')
    xfit = linspace(-2*max(abs(xmax),abs(xmin)), 2*max(abs(xmax),abs(xmin)), 1000);
    linestyle = '-';
    linewidth = 1;
elseif strcmp(fit_func, '1/x2')
    xfit = linspace(0.01, 0.2, 1000);
    linestyle = '--';
    linewidth = 0.5;
elseif strcmp(fit_func, 'critical_DTC')
    xfit = linspace(xmin, 0.2, 1000);
    linestyle = '-';
    linewidth = 0.5;
elseif strcmp(fit_func, '1/xsqrt')
    xfit = linspace(xmin*0.1, xmax*2, 1000);
    linestyle = '-';
    linewidth = 1;
elseif strcmp(fit_func, 'linear0')
    xfit = linspace(0, xmax*2, 1000);
    linestyle = '-';
    linewidth = 1;
elseif strcmp(fit_func, 'linear_s')
    xfit = linspace(-1, 1, 1000);
    linestyle = '-';
    linewidth = 1;
else
    xfit = linspace(xmin, xmax, 1000);
    linestyle = '-';
    linewidth = 2;
end


if exist('clr','var')
    if ~strcmp(clr, 'noplot')
        plot(xfit, FitRes(xfit),'linestyle', linestyle, 'linewidth', linewidth,'color', clr, 'handlevisibility', 'off')
    end
elseif strcmp(fit_func, 'DTCdecay') || strcmp(fit_func, 'logstrexp')
    plot(xfit, exp(FitRes(xfit)), 'linestyle', linestyle, 'linewidth', linewidth, 'handlevisibility', 'off')
else
    plot(xfit, FitRes(xfit), 'linestyle', linestyle, 'linewidth', linewidth, 'handlevisibility', 'off')
end


outparam = FitRes;
